docs/source/nb_examples/Continuous-Time Signals and Systems using sigsys.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pylab inline\n",
    "import sk_dsp_comm.sigsys as ss\n",
    "import scipy.signal as signal\n",
    "from IPython.display import Image, SVG"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pylab.rcParams['savefig.dpi'] = 100 # default 72\n",
    "%config InlineBackend.figure_formats=['svg'] # SVG inline viewing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to Python and the Jupyter Notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = arange(-4,4,.01)\n",
    "x = cos(2*pi*t)\n",
    "plot(t,x)\n",
    "\n",
    "grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rectangle and Triangle Pulses Defined\n",
    "Before showing more examples, consider some familiar signal primitives in your signals and systems background."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see these defined in the text see in particular Appendix F.5 (p.727) in the table of Fourier transform pairs.\n",
    "\n",
    "**Rectangle**\n",
    "\\begin{align}\n",
    "\t\\Pi\\Big(\\frac{t}{\\tau}\\Big) &= \\begin{cases}\n",
    "\t\t1, & |t| \\leq \\tau/2 \\\\\n",
    "\t\t0, & \\text{otherwise}\n",
    "\t\\end{cases}\n",
    "\\end{align}\n",
    "**Triangle**\n",
    "\\begin{align}\n",
    "\t\\Lambda\\Big(\\frac{t}{\\tau}\\Big) &= \\begin{cases}\n",
    "\t\t1-|t/\\tau|, & |t|\\leq \\tau \\\\\n",
    "\t\t0, & \\text{otherwise}\n",
    "\t\\end{cases}\n",
    "\\end{align}\n",
    "To more readily play with these function represent them numerically in Python. The module `ss.py` has some waveform primitives to help."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = arange(-5,5,.01)\n",
    "x_rect = ss.rect(t-3,2)\n",
    "x_tri = ss.tri(t+2,1.5)\n",
    "subplot(211)\n",
    "plot(t,x_rect)\n",
    "grid()\n",
    "ylabel(r'$\\Pi((t-3)/2)$');\n",
    "subplot(212)\n",
    "plot(t,x_tri)\n",
    "grid()\n",
    "xlabel(r'Time (s)')\n",
    "ylabel(r'$\\Lambda((t+2)/1.5)$');\n",
    "tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Consider an interactive version of the above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Make an interactive version of the above\n",
    "from ipywidgets import interact, interactive\n",
    "\n",
    "def pulses_plot(D1,D2,W1,W2):\n",
    "    t = arange(-5,5,.01)\n",
    "    x_rect = ss.rect(t-D1,W1)\n",
    "    x_tri = ss.tri(t-D2,W2)\n",
    "    subplot(211)\n",
    "    plot(t,x_rect)\n",
    "    grid()\n",
    "    ylabel(r'$\\Pi((t-3)/2)$');\n",
    "    subplot(212)\n",
    "    plot(t,x_tri)\n",
    "    grid()\n",
    "    xlabel(r'Time (s)')\n",
    "    ylabel(r'$\\Lambda((t+2)/1.5)$');\n",
    "    tight_layout()\n",
    "    \n",
    "interactive_plot = interactive(pulses_plot,D1 = (-3,3,.5), D2 = (-3,3,.5), W1 = (0.5,2,.25), W2 = (0.5,2,.25));\n",
    "output = interactive_plot.children[-1]\n",
    "output.layout.height = '350px'\n",
    "interactive_plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## More Signal Plotting\n",
    "The basic pulse shapes (primitives) defined in the module `ssd.py` are very useful for working Text 2.13a &d, but there are also times when you need a custom piecewise function.\n",
    "\n",
    "### Simple Cases:\n",
    "Consider plotting\n",
    "\n",
    "* $x_1(t) = \\sin(2\\pi\\cdot 5t) \\Pi((t-2)/2)$ for $0\\leq t \\leq 10$\n",
    "* $x_2(t) = \\sum_{n=-\\infty}^\\infty = \\Pi((t-5n)/1)$ for $-10 \\leq t \\leq 10$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t1 = arange(0,10+.01,.01) # arange stops one step size less than the upper limit\n",
    "x1 = sin(2*pi*5*t1)* ss.rect(t1-2,2)\n",
    "subplot(211)\n",
    "plot(t1,x1)\n",
    "xlabel(r'Time (s)')\n",
    "ylabel(r'$x_1(t)$')\n",
    "grid()\n",
    "t2 = arange(-10,10,.01)\n",
    "# Tweak mod() to take on negative values\n",
    "x2 = ss.rect(mod(t2+2.5,5)-2.5,1)\n",
    "subplot(212)\n",
    "plot(t2,x2)\n",
    "xlabel(r'Time (s)')\n",
    "ylabel(r'$x_2(t)$')\n",
    "grid()\n",
    "tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Custom Piecewise:\n",
    "A custom piecewise function is a direct and to the point way of getting a more complex function plotted. Consider plotting:\n",
    "\\begin{align}\n",
    "    x_3(t) = \\begin{cases}\n",
    "        1 + t^2, & 0\\leq t \\leq 3 \\\\\n",
    "        \\cos(2\\pi\\cdot5\\cdot t) & 3 < t \\leq 5 \\\\\n",
    "        0, & \\text{otherwise}\n",
    "    \\end{cases}\n",
    "\\end{align}\n",
    "for $-2\\leq t \\leq 6$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def x3_func(t):\n",
    "    \"\"\"\n",
    "    Create a piecewise function for plotting x3\n",
    "    \"\"\"\n",
    "    x3 = zeros_like(t)\n",
    "    for k,tk in enumerate(t):\n",
    "        if tk >= 0 and tk <= 3:\n",
    "            x3[k] = 1 + tk**2\n",
    "        elif tk > 3 and tk <= 5:\n",
    "            x3[k] = cos(2*pi*5*tk)\n",
    "    return x3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t3 = arange(-2,6+.01,.01)\n",
    "x3 = x3_func(t3)\n",
    "plot(t3,x3)\n",
    "xlabel(r'Time (s)')\n",
    "ylabel(r'$x_3(t)$')\n",
    "xlim([-2,6])\n",
    "grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "26/2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Energy and Power Signals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The general definitions are:\n",
    "\\begin{align}\n",
    "\t\tE &\\overset{\\Delta}{=} \\lim_{T\\rightarrow\\infty} \\int_{-T}^T |x(t)|^2\\, dt = \\int_{-\\infty}^\\infty |x(t)|^2\\, dt \\\\\n",
    "\t\tP &\\overset{\\Delta}{=} \\lim_{T\\rightarrow\\infty}\\frac{1}{2T} \\int_{-T}^T |x(t)|^2\\, dt\n",
    "\\end{align}\n",
    "For the case of a periodic signal, you can take the definition of $P$ above and reduce the calculation down to\n",
    "\\begin{align}\n",
    "    P = \\frac{1}{T} \\int_{t_0}^{t_0+T} |x(t)|^2\\, dt\n",
    "\\end{align}\n",
    "where $t_0$ can be any convenient value."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the waveform of Text problem 2.14b\n",
    "\\begin{align}\n",
    "    x_2(t) = \\sum_{n=-\\infty}^\\infty \\Lambda\\Big(\\frac{t-3n}{2}\\Big)\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can create an approximation to the waveform over a finite number of periods by doing a little programming:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def periodic_tri(t,tau,T,N):\n",
    "    \"\"\"\n",
    "    Approximate x2(t) by running the sum index from -N to +N.\n",
    "    The period is set by T and tau is the tri pulse width\n",
    "    parameter (base width is 2*tau).\n",
    "    \n",
    "    Mark Wickert January 2015\n",
    "    \"\"\"\n",
    "    x = zeros_like(t)\n",
    "    for n in arange(-N,N+1):\n",
    "        x += ss.tri(t-T*n,tau)\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = arange(-10,10,.001)\n",
    "x = periodic_tri(t,2,6,10)\n",
    "plot(t,x)\n",
    "plot(t,abs(x)**2)\n",
    "grid()\n",
    "#xlim([-5,5])\n",
    "xlabel(r'Time (s)')\n",
    "ylabel(r'$x_2(t)$ and $x_2^2(t)$');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the power calculation create a time array that runs over exactly one period. Below is the case for the original problem statement."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T0 = 6\n",
    "tp = arange(-T0/2,T0/2+.001,.001)\n",
    "xp = periodic_tri(tp,2,T0,5)\n",
    "plot(tp,xp)\n",
    "plot(tp,abs(xp)**2)\n",
    "legend((r'$x(t)$', r'$|x(t)|^2$'),loc='best',shadow=True)\n",
    "grid();\n",
    "xlim([-T0/2,T0/2])\n",
    "xlabel(r'Time (s)')\n",
    "ylabel(r'$x_2(t)$ and $x_2^2(t)$');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A simple numerical approximation to the integral\n",
    "\\begin{align}\n",
    "    P = \\frac{1}{T}\\int_0^T |x_b(t)|^2\\, dt\n",
    "\\end{align}\n",
    "is shown below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Power calculation\n",
    "Px2 = (1/T0)*sum(xp**2)*.001 # rectangular partitions for integral\n",
    "print('Power estimate via numerical integration: %2.4f W' % Px2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Power in the Sum of Two Sinusoids"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The problem is what is the power in the signal\n",
    "\\begin{align}\n",
    "    x(t) = A_1 \\cos(\\omega_1 t +\\phi_1) + A_2 \\cos(\\omega_2 t + \\phi_2),\\ -\\infty < t < \\infty\n",
    "\\end{align}\n",
    "Since we are not certain that $x(t)$ is periodic, the power calculation requires that we form\n",
    "\\begin{align}\n",
    "    P_x = \\lim_{T\\rightarrow\\infty} \\frac{1}{T} \\int_{-T/2}^{T/2} |x(t)|^2\\, dt = \\langle |x(t)|^2\\rangle\n",
    "\\end{align}\n",
    "* Rather that just jumping in and making a mess, consider first the expansion of $|x(t)|^2 = x^2(t)$:\n",
    "\\begin{align}\n",
    "    x^2(t) &= \\frac{A_1^2}{2}\\big[1+\\cos(2\\omega_1 t + \\phi_1)\\big] + \\frac{A_2^2}{2}\\big[1+\\cos(2\\omega_2 t + \\phi_2)\\big] \\\\\n",
    "    &\\quad + 2\\frac{A_1 A_2}{2}\\Big\\{\\cos[(\\omega_1 + \\omega_2)t + (\\phi_1+\\phi_2)\\big] + \\cos[(\\omega_1 - \\omega_2)t + (\\phi_1-\\phi_2)\\big]\\Big\\}\n",
    "\\end{align}\n",
    "* The time average operator is linear, so we consider $\\langle\\ \\ \\rangle$ operating on each term of the above independently\n",
    "* For $\\omega_1 \\neq \\omega_2$, the first two terms yield $A_1^2/2$ and $A_2^2/2$ respectively\n",
    "* The last term requires some thinking, but as long as $\\omega_1 \\neq \\omega_2$ the times average of $\\cos[(\\omega_1 + \\omega_2)t + (\\phi_1+\\phi_2)]$ and $\\cos[(\\omega_1 - \\omega_2)t + (\\phi_1-\\phi_2)$], the two terms respectively are each zero!\n",
    "* Finally,\n",
    "\\begin{align}\n",
    "    P_x = \\frac{A_1^2}{2} + \\frac{A_2^2}{2}\n",
    "\\end{align}\n",
    "* When the frequencies are equal, then you can combine the terms using trig identities (recall the phasor addition formula from ECE 2610\n",
    "\\begin{align}\n",
    "    x(t) = A\\cos(\\omega t + \\phi)\n",
    "\\end{align}\n",
    "where $\\omega = \\omega_1 = \\omega_2$ and\n",
    "\\begin{align}\n",
    "    Ae^{j\\phi} = A_1e^{j\\phi_1} + A_2 e^{j\\phi_2}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = arange(-10,10,.001)\n",
    "x1 = 4*cos(2*pi*10*t)\n",
    "x2 = 3*cos(2*pi*3.45*t+pi/9)\n",
    "plot(t,x1)\n",
    "plot(t,x2)\n",
    "plot(t,x1+x2)\n",
    "grid()\n",
    "xlabel(r'Time (s)')\n",
    "ylabel(r'Amplitude')\n",
    "legend((r'$x_1(t)$', r'$x_2(t)$', r'$x_1(t)+x_2(t)$'),loc='best',shadow=True)\n",
    "xlim([-.1,.1]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Power calculations: %3.2f, %3.2f, %3.2f' \\\n",
    "      % (var(x1),var(x2),var(x1+x2)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Theory: %3.2f, %3.2f, %3.2f' \\\n",
    "      % (4**2/2,3**2/2,4**2/2+3**2/2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fourier Series and Line Spectra Plotting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Being able to easily plot the line spectra of periodic signals will hopefully enhance your understanding. The module `ss.py` contains the function `ss.line_spectra()` for this purpose. The function assumes that the Fourier coefficients, $X_n$ are available for a real signal $x(t)$. The function plots line spectra as:\n",
    "* The two-sided magnitude spectra\n",
    "* The two-sided magnitude spectra in dB with an adjustable floor level in dB\n",
    "* The two-sided phase spectra in radians\n",
    "* The one-sided line spectra corresponding to the three cases listed immediately above\n",
    "Examples are given below for the case of a simple pulse train and then for a trapezoidal pulse train. IN the case of the trapezoidal pulse train the underlying Fourier coefficients are obtained numerically using the FFT as described in the course notes.\n",
    "\n",
    "A fundamental requirement in using `ss.line_spectra()` is to beable to supply the coefficients starting with the DC term coefficient $X_0$ and moving up to the $N$th harmonic. Before plotting the pulse train line spectra I first describe a *helper* function for visualizing the pulse train waveform."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pulse Train"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pulse_train(Np,fs,tau,t0):\n",
    "    \"\"\"\n",
    "    Generate a discrete-time approximation to a continuous-time\n",
    "    pulse train signal. Amplitude values are [0,1]. Scale and offset\n",
    "    later if needed.\n",
    "    \n",
    "    Inputs\n",
    "    ------\n",
    "     Np = number of periods to generate\n",
    "     fs = samples per period\n",
    "    tau = duty cycle\n",
    "     t0 = pulse delay time relative to first rising edge at t = 0\n",
    "    \n",
    "    Return\n",
    "    ------\n",
    "    t = time axis array\n",
    "    x = waveform\n",
    "    \n",
    "    Mark Wickert, January 2015\n",
    "    \"\"\"\n",
    "    t = arange(0,Np*fs+1,1)/fs #time is normalized to make period T0 = 1.0\n",
    "    x = zeros_like(t)\n",
    "    # Using a brute force approach, just fill x with the sample values\n",
    "    for k,tk in enumerate(t):\n",
    "        if mod(tk-t0,1) <= tau and mod(tk-t0,1) >= 0:\n",
    "            x[k] = 1\n",
    "    return t,x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tau = 1/8; fs = 8*16; t0 = 0 # note t0 = tau/2\n",
    "subplot(211)\n",
    "t,x = pulse_train(4,fs,tau,t0)\n",
    "plot(t,x) # Just a plot of xa(t)\n",
    "ylim([-.1,1.1])\n",
    "grid()\n",
    "ylabel(r'$x_a(t)$')\n",
    "title(r'Pulse Train Signal: (top) $x_a(t)$, (bot) $x_b(t) = 1-x_a(t)$');\n",
    "subplot(212)\n",
    "t,x = pulse_train(4,fs,tau,t0)\n",
    "plot(t,1-x) # Note here y(t) = 1 - x(t), a special case of \n",
    "ylim([-.1,1.1]) # y(t) = A + B*x(t) in the notes\n",
    "grid()\n",
    "xlabel(r'Time ($t/T_0$)')\n",
    "ylabel(r'$x_b(t)$');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: Pulse Train Line Spectra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the case of pulse train having the initial pulse starting at $t=0$, i.e.,\n",
    "\\begin{align}\n",
    "    x(t) = \\sum_{k=-\\infty}^\\infty A\\cdot \\Pi\\left(\\frac{t-\\tau/2-kT_0}{\\tau}\\right),\n",
    "\\end{align}\n",
    "the Fourier coefficient are given by\n",
    "\\begin{align}\n",
    "    X_n = A\\cdot\\frac{\\tau}{T_0}\\cdot\\text{sinc}(nf_0\\tau)\\cdot\\exp(-j2\\pi n f_0t_0)\n",
    "\\end{align}\n",
    "where $f_0 = 1/T_0$ is the fundamental frequency and here $t_0 = \\tau/2$.\n",
    "\n",
    "Line spectra plotting is shown below for this case. If the pulse train should be shifted in time to some other orientation, then the phase plot will change, as the included $\\exp(j2\\pi n f_0 t_0)$ term will be different.\n",
    "\n",
    "**Note:** The pulse train function define above is slightly different from the pulse train defined in the book and shown in mathematical form as $x(t)$ just above in this cell. The function `pulse_train()` has the first pulse starting exactly at $t=0$. To move the pule train right or left on the time axis, you can use the function parameter `t0`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = arange(0,25+1) # Get 0 through 25 harmonics\n",
    "tau = 0.125; f0 = 1; A = 1;\n",
    "Xn = A*tau*f0*sinc(n*f0*tau)*exp(-1j*2*pi*n*f0*tau/2)\n",
    "# Xn = -Xn # Convert the coefficients from xa(t) t0 xb(t)\n",
    "# Xn[0] += 1\n",
    "figure(figsize=(6,2))\n",
    "f = n # Assume a fundamental frequency of 1 Hz so f = n\n",
    "ss.line_spectra(f,Xn,mode='mag',sides=2,fsize=(6,2))\n",
    "xlim([-25,25]);\n",
    "#ylim([-50,10])\n",
    "figure(figsize=(6,2))\n",
    "ss.line_spectra(f,Xn,mode='phase',fsize=(6,2))\n",
    "xlim([-25,25]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: Trapezoidal Pulse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the line spectra of a finite rise and fall time pulse train is of practical interest. The function `trap_pulse()` allows you first visualize one period of the trapezoidal pulse train, and then use this waveform in obtaining numerically the Fourier coefficients of this signal. PLotting the corresponding line spectra follows.\n",
    "\n",
    "A point to be main is that by slowing down the edges (rise time/fall time) of the pulse train the amplitude of the harmonics falls off more rapidly. When considering the clock speed in todays PCs this can be a good thing as harmonic emission is an issue."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def trap_pulse(N,tau,tr):\n",
    "    \"\"\"\n",
    "    xp = trap_pulse(N,tau,tr)\n",
    "    \n",
    "    Mark Wickert, January 2015\n",
    "    \"\"\"\n",
    "    n = arange(0,N)\n",
    "    t = n/N\n",
    "    xp = zeros(len(t))\n",
    "    # Assume tr and tf are equal\n",
    "    T1 = tau + tr\n",
    "    # Create one period of the trapezoidal pulse waveform\n",
    "    for k in n:\n",
    "        if t[k] <= tr:\n",
    "            xp[k] = t[k]/tr\n",
    "        elif (t[k] > tr and t[k] <= tau):\n",
    "            xp[k] = 1\n",
    "        elif (t[k] > tau and t[k] < T1):\n",
    "            xp[k] = -t[k]/tr + 1 + tau/tr;\n",
    "        else:\n",
    "            xp[k] = 0\n",
    "    return xp, t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let $\\tau = 1/8$ and $t_r = 1/20$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# tau = 1/8, tr = 1/20\n",
    "N = 1024\n",
    "xp,t = trap_pulse(N,1/8,1/20)\n",
    "Xp = fft.fft(xp)\n",
    "figure(figsize=(6,2))\n",
    "plot(t,xp)\n",
    "grid()\n",
    "title(r'Spectra of Finite Risetime Pulse Train: $\\tau = 1/8$ $t_r = 1/20$')\n",
    "ylabel(r'$x(t)$')\n",
    "xlabel('Time (s)')\n",
    "f = arange(0,N/2)\n",
    "ss.line_spectra(f[0:25],Xp[0:25]/N,'magdB',floor_dB=-80,fsize=(6,2))\n",
    "ylabel(r'$|X_n| = |X(f_n)|$ (dB)');\n",
    "#% tau = 1/8, tr = 1/10\n",
    "xp,t = trap_pulse(N,1/8,1/10)\n",
    "Xp = fft.fft(xp)\n",
    "figure(figsize=(6,2))\n",
    "plot(t,xp)\n",
    "grid()\n",
    "title(r'Spectra of Finite Risetime Pulse Train: $\\tau = 1/8$ $t_r = 1/10$')\n",
    "ylabel(r'$x(t)$')\n",
    "xlabel('Time (s)')\n",
    "ss.line_spectra(f[0:25],Xp[0:25]/N,'magdB',floor_dB=-80,fsize=(6,2))\n",
    "ylabel(r'$|X_n| = |X(f_n)|$ (dB)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the edge speed slowed down it is clear that the harmonics drop off faster."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fourier Transforms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Fourier transfrom definition is:\n",
    "\\begin{align}\n",
    "    X(f) &= \\int_{-\\infty}^\\infty x(t)\\ e^{-j2\\pi ft}\\, dt \\\\\n",
    "    x(t) &= \\int_{-\\infty}^\\infty X(f)\\, e^{j2\\pi ft}\\, df\n",
    "\\end{align}\n",
    "\n",
    "A numerical approximation to the Fourier transform is possible using the FFT, or more conveniently using the function `freqz()` from the package `scipy.signal`. A helper function to abstract some of the digital signal processing details is `f, X = FT_approx(x,dt,Nfft)`. The function is now part of `sigsys.py` with name change to `ft_approx()`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def FT_approx(x,t,Nfft):\n",
    "    '''\n",
    "    Approximate the Fourier transform of a finite duration\n",
    "    signal using scipy.signal.freqz()\n",
    "    \n",
    "    Inputs\n",
    "    ------\n",
    "       x = input signal array\n",
    "       t = time array used to create x(t)\n",
    "    Nfft = the number of frdquency domain points used to \n",
    "           approximate X(f) on the interval [fs/2,fs/2], where\n",
    "           fs = 1/Dt. Dt being the time spacing in array t\n",
    "    \n",
    "    Return\n",
    "    ------\n",
    "    f = frequency axis array in Hz\n",
    "    X = the Fourier transform approximation (complex)\n",
    "    \n",
    "    Mark Wickert, January 2015\n",
    "    '''\n",
    "    fs = 1/(t[1] - t[0])\n",
    "    t0 = (t[-1]+t[0])/2 # time delay at center\n",
    "    N0 = len(t)/2 # FFT center in samples\n",
    "    f = arange(-1/2,1/2,1/Nfft)\n",
    "    w, X = signal.freqz(x,1,2*pi*f)\n",
    "    X /= fs # account for dt = 1/fs in integral\n",
    "    X *= exp(-1j*2*pi*f*fs*t0)# time interval correction\n",
    "    X *= exp(1j*2*pi*f*N0)# FFT time interval is [0,Nfft-1]\n",
    "    F = f*fs\n",
    "    return F, X"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: Rectangular Pulse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As as simple starting point example, consider $x(t) = \\Pi(t\\tau)$. The well known result for the Fourier transfrom (FT) is:\n",
    "\\begin{align}\n",
    "    X(f) = \\mathcal{F}\\left\\{\\Pi\\left(\\frac{t}{\\tau}\\right)\\right\\} = \\tau\\,\\text{sinc}(f\\tau)\n",
    "\\end{align}\n",
    "We now use the above defined `FT_approx()` to obtain a numerical approximation to the FT of the rectangular pulse.\n",
    "\n",
    "**Tips:**\n",
    "* Make sure the signal is well contained on the time interval used to generate $x(t)$\n",
    "* Make sure the sampling rate, one over the sample spacing, is adequate to represent the signal spectrum\n",
    "* From sampling theory, the reange of frequencies represented by the spectrum estimate will be $f_s/2 \\leq f < f_s/2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fs = 100 # sampling rate in Hz\n",
    "tau = 1\n",
    "t = arange(-5,5,1/fs)\n",
    "x0 = ss.rect(t-.5,tau)\n",
    "figure(figsize=(6,5))\n",
    "subplot(311)\n",
    "plot(t,x0)\n",
    "grid()\n",
    "ylim([-0.1,1.1])\n",
    "xlim([-2,2])\n",
    "title(r'Exact Waveform')\n",
    "xlabel(r'Time (s)')\n",
    "ylabel(r'$x_0(t)$');\n",
    "\n",
    "# FT Exact Plot\n",
    "fe = arange(-10,10,.01)\n",
    "X0e = tau*sinc(fe*tau)\n",
    "subplot(312)\n",
    "plot(fe,abs(X0e))\n",
    "#plot(f,angle(X0))\n",
    "grid()\n",
    "xlim([-10,10])\n",
    "title(r'Exact Spectrum Magnitude')\n",
    "xlabel(r'Frequency (Hz)')\n",
    "ylabel(r'$|X_0e(f)|$');\n",
    "\n",
    "# FT Approximation Plot\n",
    "f,X0 = ss.ft_approx(x0,t,4096)\n",
    "subplot(313)\n",
    "plot(f,abs(X0))\n",
    "#plot(f,angle(X0))\n",
    "grid()\n",
    "xlim([-10,10])\n",
    "title(r'Approximation Spectrum Magnitude')\n",
    "xlabel(r'Frequency (Hz)')\n",
    "ylabel(r'$|X_0(f)|$');\n",
    "tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: Text Problem 2.31a Drill Down"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this problem you are given\n",
    "\\begin{align}\n",
    "    x_1(t) = \\Pi\\left(\\frac{t+1/2}{1}\\right) - \\Pi\\left(\\frac{t-1/2}{1}\\right)\n",
    "\\end{align}\n",
    "The Fourier transfrom of this signal can be found using *linearity* and the *time delay* theorems.\n",
    "\\begin{align}\n",
    "    X_1(f) &= \\mathcal{F}\\left\\{\\Pi\\left(\\frac{t+1/2}{1}\\right) - \\Pi\\left(\\frac{t-1/2}{1}\\right)\\right\\} \\\\\n",
    "    &= \\text{sinc}(f)\\cdot\\left[e^{j2\\pi f\\cdot 1/2} - e^{-j2\\pi f\\cdot 1/2}\\right]\\times\\frac{2j}{2j} \\\\\n",
    "    &= 2j\\ \\text{sinc}(f)\\cdot\\sin(\\pi f)\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fs = 100\n",
    "t = arange(-5,5,1/fs)\n",
    "x1 = ss.rect(t+1/2,1)-ss.rect(t-1/2,1)\n",
    "subplot(211)\n",
    "plot(t,x1)\n",
    "grid()\n",
    "ylim([-1.1,1.1])\n",
    "xlim([-2,2])\n",
    "xlabel(r'Time (s)')\n",
    "ylabel(r'$x_1(t)$');\n",
    "fe = arange(-10,10,.01)\n",
    "X1e = 2*1j*sinc(fe)*sin(pi*fe)\n",
    "f,X1 = ss.ft_approx(x1,t,4096)\n",
    "subplot(212)\n",
    "plot(f,abs(X1))\n",
    "plot(fe,abs(X1e))\n",
    "#plot(f,angle(X1))\n",
    "legend((r'Num Approx',r'Exact'),loc='best')\n",
    "grid()\n",
    "xlim([-10,10])\n",
    "xlabel(r'Frequency (Hz)')\n",
    "ylabel(r'$|X_1(f)|$');\n",
    "tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Notice the numerical approximation and exact spectral plots overlay one another"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: Modulation Theorem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the modulation theorem, which is extremely important to communications theory:\n",
    "\\begin{align}\n",
    "    y(t) &= x(t)\\cdot\\cos(2\\pi f_0 t) \\\\\n",
    "    Y(f) &= \\frac{1}{2}\\left[X(f-f_0) + X(f+f_0)\\right]\n",
    "\\end{align}\n",
    "Here we will use a triangle pulse for $x(t)$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fs = 100 # sampling rate in Hz\n",
    "tau = 1\n",
    "t = arange(-5,5,1/fs)\n",
    "x3 = ss.tri(t,tau)\n",
    "y = x3*cos(2*pi*10*t)\n",
    "subplot(211)\n",
    "plot(t,x3)\n",
    "plot(t,y)\n",
    "grid()\n",
    "ylim([-1.1,1.1])\n",
    "xlim([-2,2])\n",
    "legend((r'$x_3(t)$', r'$y(t)$'),loc='lower right',shadow=True)\n",
    "title(r'Time Domain: $x_3(t)$ and $y(t)=x_3(t)\\cos(2\\pi\\cdot 5\\cdot t)$')\n",
    "xlabel(r'Time (s)')\n",
    "ylabel(r'$y(t)$');\n",
    "f,Y = ss.ft_approx(y,t,4096)\n",
    "subplot(212)\n",
    "plot(f,abs(Y))\n",
    "#plot(f,angle(X0))\n",
    "grid()\n",
    "title(r'Frequency Domain: $Y(f)$')\n",
    "xlim([-15,15])\n",
    "xlabel(r'Frequency (Hz)')\n",
    "ylabel(r'$|Y(f)|$');\n",
    "tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: Representing a Bandlimited Signal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We know that in theory a bandlimited signal can only be generated from a signal having infinite duration. Specifically, a signal with rectangular spectrum has Fourier transfrom pair:\n",
    "\\begin{align}\n",
    "    x(t) = 2W\\text{sinc}(2Wt) \\overset{\\mathcal{F}}{\\Leftrightarrow} \\Pi\\left(\\frac{f}{2W}\\right) = X(f)\n",
    "\\end{align}\n",
    "In a simulation we expect to have troubles modeling the finite duration aspects of the signal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fs = 100 # sampling rate in Hz\n",
    "W = 5\n",
    "t = arange(-5,5,1/fs)\n",
    "x4 = 2*W*sinc(2*W*t)\n",
    "figure(figsize=(6,2))\n",
    "plot(t,x4)\n",
    "grid()\n",
    "#ylim([-1.1,1.1])\n",
    "xlim([-2,2])\n",
    "title(r'Time Domain: $x_4(t),\\ W = 5$ Hz')\n",
    "xlabel(r'Time (s)')\n",
    "ylabel(r'$x_4(t)$');\n",
    "f,X4 = ss.ft_approx(x4,t,4096)\n",
    "figure(figsize=(6,2))\n",
    "plot(f,abs(X4))\n",
    "grid()\n",
    "title(r'Frequency Domain: $X_4(f)$')\n",
    "xlim([-10,10])\n",
    "xlabel(r'Frequency (Hz)')\n",
    "ylabel(r'$|X_4(f)|$');\n",
    "figure(figsize=(6,2))\n",
    "plot(f,20*log10(abs(X4)))\n",
    "grid()\n",
    "title(r'Frequency Domain: $X_4(f)$ in dB')\n",
    "ylim([-50,5])\n",
    "xlim([-10,10])\n",
    "xlabel(r'Frequency (Hz)')\n",
    "ylabel(r'$|X_4(f)|$ (dB)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note:** The dB version (last plot) reveals that the first sidelobes of the spectrum are only down ~21dB. Increasing the length of the time window will not help. The spectral side lobes will become more tightly packed, but the first sidelobe will still be down only 21dB. With other pulse shapes in the time domain, i.e., note simply a truncted $\\text{sinc}()$ function reduced sidelobes can be obtained."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Convolution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* The convolution of two signals $x_1(t)$ and $x_2(t)$ is defined as\n",
    "\\begin{align}\n",
    "\tx(t) &= x_1(t)\\ast x_2(t) = \\int_{-\\infty}^\\infty x_1(\\lambda)x_2(t-\\lambda)\\, d\\lambda \\\\\n",
    "\t&= x_2(t)\\ast x_1(t) = \\int_{-\\infty}^\\infty x_2(\\lambda)x_1(t-\\lambda)\\, d\\lambda\n",
    "\\end{align}\n",
    "* A special convolution case is $\\delta(t-t_0)$\n",
    "\\begin{align}\n",
    "\t\t\\delta(t-t_0)\\ast x(t) &= \\int_{-\\infty}^\\infty \\delta(\\lambda-t_0)x(t-\\lambda)\\, d\\lambda \\\\\n",
    "\t\t&= x(t-\\lambda)\\big|_{\\lambda=t_0} = x(t-t_0)\n",
    "\\end{align}\n",
    "You can experiment with the convolution integral numerically using `ssd.conv_integral()` found in the module `ssd.py`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = arange(-2,2.001,.001)\n",
    "p1 = ss.rect(t,1)\n",
    "p2 = ss.rect(t,3)\n",
    "y,ty = ss.conv_integral(p1,t,p2,t)\n",
    "plot(ty,y)\n",
    "ylim([-.01,1.01])\n",
    "grid()\n",
    "xlabel(r'Time (s)')\n",
    "ylabel(r'$y(t)$');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For convolutions involving semi-infinite signals, such as $u(t)$, you can tell `ssd.conv_integral()` about this via the optional extent argument. See the function help using\n",
    "```python\n",
    "ss.conv_integral?\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Consider a pulse convolved with an exponential ('r' type extent)\n",
    "tx = arange(-1,8,.01)\n",
    "x = ss.rect(tx-2,4) # pulse starts at t = 0\n",
    "h = 4*exp(-4*tx)*ss.step(tx)\n",
    "y,ty = ss.conv_integral(x,tx,h,tx,extent=('f','r')) # note extents set\n",
    "plot(ty,y) # expect a pulse charge and discharge waveform\n",
    "grid()\n",
    "title(r'$\\Pi((t-2)/4)\\ast 4 e^{-4t} u(t)$')\n",
    "xlabel(r'Time (s)')\n",
    "ylabel(r'$y(t)$');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Spectrum of PN Sequence (exact)\n",
    "The cell below is a copy of the earlier pulse train line spectra example. Use this as a template to create the solution to the PN code problem of HW 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = arange(0,25+1) # Get 0 through 25 harmonics\n",
    "tau = 0.125; f0 = 1; A = 1;\n",
    "Xn = A*tau*f0*sinc(n*f0*tau)*exp(-1j*2*pi*n*f0*tau/2)\n",
    "# Xn = -Xn # Convert the coefficients from xa(t) t0 xb(t)\n",
    "# Xn[0] += 1\n",
    "figure(figsize=(6,2))\n",
    "f = n # Assume a fundamental frequency of 1 Hz so f = n\n",
    "ss.line_spectra(f,Xn,mode='mag',sides=2,fsize=(6,2))\n",
    "xlim([-25,25]);\n",
    "#ylim([-50,10])\n",
    "figure(figsize=(6,2))\n",
    "ss.line_spectra(f,Xn,mode='phase',fsize=(6,2))\n",
    "xlim([-25,25]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Spectrum of PN Sequence (approx)\n",
    "The code below approximates the PSD of the PN code using a numerical approximation to the Fourier coefficients, $X_n$. This development may be useful for the lab, as you can esily change the waveform level without having to rework the theory.\n",
    "\n",
    "The approach taken here to create one period of the PN waveform at 10 samples per bit. The line containing the function `ss.upsample()` converts the bit sequence into a waveform by upsampling and filtering with a rectangular pulse shape (`ones(10)`). The function `ss.fs_coeff()` numerically calculates the $X_n$'s. To plot the PSD from the Fourier coefficients we use\n",
    "\n",
    "$$\n",
    "    S_x(f) = \\sum_{n=-\\infty}^\\infty |X_n|^2 \\delta(f-nf_0)\n",
    "$$\n",
    "\n",
    "where $f_0$ in this case is $1/(MT_0)$ with $T_0$ beging the bit period and $M$ the code period in bits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_PN4 = ss.m_seq(4)\n",
    "x = signal.lfilter(ones(10),1,ss.upsample(x_PN4,10))\n",
    "t = arange(0,len(x))/10\n",
    "figure(figsize=(6,2))\n",
    "plot(t,x);\n",
    "title(r'Time Domain and PSD of $M=15$ PN Code with $T = 1$')\n",
    "xlabel(r'Time (s)')\n",
    "ylabel(r'x(t)')\n",
    "axis([0,15,-0.1,1.1]);\n",
    "grid()\n",
    "# 10 samples/bit so 150 samples/period\n",
    "# harmonics spaced by 1/(15*T) = 1/15\n",
    "Xk,fk = ss.fs_coeff(x,45,1/15)\n",
    "ss.line_spectra(fk,Xk,'magdB',lwidth=2.0,floor_dB=-50,fsize=(6,2))\n",
    "xlim([-3,3])\n",
    "ylabel(r'$|X_n| = |X(f_n)|$ (dB)');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Line spacing\n",
    "1/15"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sk_dsp_comm.digitalcom as dc\n",
    "y_PN5_bits = ss.PN_gen(10000,5)\n",
    "# Convert to waveform level shifted to +/-1 amplitude\n",
    "y = 2*signal.lfilter(ones(10),1,ss.upsample(y_PN5_bits,10))-1\n",
    "# Find the time averaged autocorrelation function normalized\n",
    "# to have a peak amplitude of 1\n",
    "Ry,tau = dc.xcorr(y,y,400)\n",
    "# We know Ry is real so strip small imag parts from FFT-based calc\n",
    "Ry = Ry.real"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fs = 10\n",
    "t = arange(len(y))/fs\n",
    "plot(t[:500],y[:500])\n",
    "title(r'PN Waveform for 5 Stages (Period $2^5 -1 = 31$ bits)')\n",
    "ylabel(r'Amplitude')\n",
    "xlabel(r'Bits (10 samples/bit)')\n",
    "grid();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tau_s = tau/10\n",
    "figure(figsize=(6,2))\n",
    "plot(tau_s,Ry)\n",
    "title(r'Autocorrelation and PSD Estimates for $M=31$ with $T = 1$')\n",
    "xlabel(r'Autocorrelation Lag $\\tau$ (s)')\n",
    "ylabel(r'$R_y(\\tau)$')\n",
    "grid();\n",
    "figure(figsize=(6,2))\n",
    "psd(y,2**12,10)\n",
    "xlabel(r'Frequency (Hz)')\n",
    "ylabel(r'$S_y(f)$ (dB)')\n",
    "#xlim([0,.002]);\n",
    "ylim([-30,20]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Lab Tip: PN Generator Coded at the Bit Level\n",
    "In Lab 2 of ECE 4670 a C/C++ version of a PN generator is implemented to run the ARM `mbed` LPC 1768 microcontroller (https://www.sparkfun.com/products/9564). At the heart of this code is:\n",
    "```C\n",
    "// Globals defined as unsigned int\n",
    "tap1 -= 1;\n",
    "tap2 -= 1;\n",
    "mask1 = 0x1 << (tap1);\n",
    "mask2 = 0x1 << (tap2);\n",
    "bit = 0x0;\n",
    "sync = 0x0;\n",
    "\n",
    "void gen_PN() {\n",
    "    my_pin5 = bit;\n",
    "    my_pin6 = synch_bit;\n",
    "    led2 = bit;\n",
    "    led3 = synch_bit;\n",
    "    if (clk_state == 0x1)\n",
    "    {\n",
    "    // Advance m-sequence generator by one bit\n",
    "    // XOR tap1 and tap2 SR values and feedback to input\n",
    "        fb = ((sr & mask1)>> tap1) ^ ((sr & mask2) >> tap2);\n",
    "        sr = (sr << 1) + fb;\n",
    "        bit = sr & 0x1;\n",
    "        // Use random number generator in place of m-sequence bits\n",
    "        if (DIP_sw4)\n",
    "        {\n",
    "            bit = rand_int() & 0x1;\n",
    "        }\n",
    "        clk_state = 0x0;\n",
    "        // See if all 1's condition exists in SR\n",
    "        if ((sr & synch) == synch) {\n",
    "            synch_bit = 0x1;\n",
    "        }\n",
    "        else\n",
    "        {\n",
    "            synch_bit = 0x0;  \n",
    "        }\n",
    "    }\n",
    "    else\n",
    "    {\n",
    "        if (DIP_sw1) bit = !bit;\n",
    "        clk_state = 0x1;\n",
    "    }\n",
    "}\n",
    "```\n",
    "The data type is `unsigned int`, which on the mbed is `uint16_t`, an unsigned 16-bit integer. A single unsigned integer is used as a 16-bit shift register with the LSB, furthest bit to the right, used to represent the first register stage. The shift register is advanced using a left shift `<<` bitwise operation. We can code this Python almost directly, as shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class bitwise_PN(object):\n",
    "    \"\"\"\n",
    "    Implement a PN generator using bitwise manipulation for\n",
    "    the shift register. The LSB holds b0 and bits are shifted left.\n",
    "              +----+----+----+----+----+----+----+\n",
    "         sr = |bM-1| .. |bM-k| .. | b2 | b1 | b0 |\n",
    "              +----+----+----+----+----+----+----+\n",
    "                 |         |\n",
    "    Feedback:(tap1-1)   (tap2-1)    Shift left using <<\n",
    "\n",
    "    Mark Wickert February 2017\n",
    "    \"\"\"\n",
    "    def __init__(self,tap1,tap2,Nstage,sr_initialize):\n",
    "        \"\"\"\n",
    "        Initialize the PN generator object\n",
    "        \"\"\"\n",
    "        self.tap1 = tap1 - 1\n",
    "        self.tap2 = tap2 - 1\n",
    "        self.mask1 = 0x1 << (tap1 - 1) # to select bit of interest\n",
    "        self.mask2 = 0x1 << (tap2 - 1) # to select bit of interest\n",
    "        self.Nstage = Nstage\n",
    "        self.period = 2**Nstage - 1\n",
    "        self.sr = sr_initialize\n",
    "        self.bit = 0\n",
    "        self.sync_bit = 0\n",
    "    \n",
    "    def clock_PN(self):\n",
    "        '''\n",
    "        Method to advance m-sequence generator by one bit\n",
    "        XOR tap1 and tap2 SR values and feedback to input\n",
    "        '''\n",
    "        fb = ((self.sr & self.mask1)>> self.tap1) ^ \\\n",
    "            ((self.sr & self.mask2) >> self.tap2)\n",
    "        self.sr = (self.sr << 1) + fb\n",
    "        self.sr = self.sr & self.period # set MSBs > Nstage to 0\n",
    "        self.bit = self.sr & 0x1 # output LSB from SR \n",
    "        # See if all 1's condition exits in SR, if so output a synch pulse\n",
    "        if ((self.sr & self.period) == self.period):\n",
    "            self.sync_bit = 0x1\n",
    "        else:\n",
    "            self.sync_bit = 0x0  \n",
    "        print('output = %d, sr contents = %s, sync bit = %d' \\\n",
    "             % (self.bit, binary(self.sr, self.Nstage), self.sync_bit))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# A simple binary format display function which shows\n",
    "# leading zeros to a fixed bit width\n",
    "def binary(num, length=8):\n",
    "    return format(num, '#0{}b'.format(length + 2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PN1 = bitwise_PN(10,7,10,0x1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PN1.clock_PN()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sr initial condition\n",
    "sr = 0b1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Nout = 20\n",
    "x_out = zeros(Nout)\n",
    "s_out = zeros(Nout)\n",
    "PN1 = bitwise_PN(3,2,3,0x1)\n",
    "for k in range(Nout):\n",
    "    PN1.clock_PN()\n",
    "    x_out[k] = PN1.bit\n",
    "    s_out[k] = PN1.sync_bit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stem(x_out)\n",
    "stem(0.2*s_out,markerfmt = 'ro')\n",
    "ylim([0,1.1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cross Correlation and Signal Delay"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The idea of the autocorrelation function can be extended to the cross correlation, that is the correlation or likeness between two signals, say $x(t)$ and $y(t)$. Define\n",
    "\\begin{align}\n",
    "    R_{xy}(\\tau) = \\langle x(t)y(t+\\tau)\\rangle = \\lim_{T\\rightarrow\\infty} \\frac{1}{2T}\\int_{-T}^T x(t)y(t+\\tau)\\, dt\n",
    "\\end{align}\n",
    "Consider a simulation example using `dc.xcorr(x,t,lags)`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sk_dsp_comm.digitalcom as dc\n",
    "x_PN4_bits = ss.PN_gen(10000,6)\n",
    "# Convert to waveform level shifted to +/-1 amplitude\n",
    "x_s = 2*signal.lfilter(ones(10),1,ss.upsample(x_PN4_bits,10))-1\n",
    "# Form a delayed version of x_S\n",
    "T_D = 35 # 35 sample delay\n",
    "y_s = signal.lfilter(concatenate((zeros(T_D),array([1]))),1,x_s)\n",
    "figure(figsize=(6,2))\n",
    "plot(x_s[:200])\n",
    "plot(y_s[:200])\n",
    "ylim([-1.1,1.1])\n",
    "title(r'Delayed and Undelayed Signals for $T_D = 35$ Samples')\n",
    "xlabel(r'Samples (10/PN bit)')\n",
    "ylabel(r'$x_s(t)$ and $y_s(t)$')\n",
    "grid();\n",
    "# Find the time averaged autocorrelation function normalized\n",
    "# to have a peak amplitude of 1\n",
    "Ryx,tau = dc.xcorr(y_s,x_s,200) #note order change\n",
    "# We know Ryx is real\n",
    "Ryx = Ryx.real\n",
    "tau_s = tau/10\n",
    "figure(figsize=(6,2))\n",
    "plot(tau_s,Ryx)\n",
    "title(r'Cross Correlation for $M=4$ with $T = 1$ and Delay 35 Samples')\n",
    "xlabel(r'Autocorrelation Lag $\\tau$ (s)')\n",
    "ylabel(r'$R_{yx}(\\tau)$')\n",
    "grid();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Spectral Containment Bandwidth (text problem 2.55)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In text problem 2.55 you are asked to find the 90% energy contain bandwidth of a signal $x_i(t)$. Specifically you are to find the lowpass or one-sided bandwidth of a baseband signal such that 90% of the total signal energy is contained in the bandwidth, $B_{90}$. You obtain $B_{90}$ by solving the following equation\n",
    "\\begin{align}\n",
    "    0.9 = \\frac{0.9 E_\\text{total}}{E_\\text{total}} = \\frac{\\int_{-B_{90}}^{B_{90}} G(f) df}{\\int_{-\\infty}^\\infty G(f) df} = \\frac{2\\int_0^{B_{90}} G(f) df}{2\\int_0^\\infty G(f) df} = \\frac{\\int_0^{B_{90}} G(f) df}{\\int_0^\\infty G(f) df},\n",
    "\\end{align}\n",
    "where $G(f) = |X_i(f)|^2$ is the energy spectral density of $x_i(t)$.\n",
    "\n",
    "For parts (c) and (d) the problem states you need to perform numerical integration."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In an exalier example found in this notebook I found the Fourier transform of\n",
    "\\begin{align}\n",
    "    x(t) = \\Pi\\left(\\frac{t-\\tau/2}{\\tau}\\right) - \\Pi\\left(\\frac{t+\\tau/2}{\\tau}\\right)\n",
    "\\end{align}\n",
    "to be\n",
    "\\begin{align}\n",
    "    X(f) &= 2j\\ \\text{sinc}(f\\tau)\\cdot\\sin(\\pi f\\tau)\n",
    "\\end{align}\n",
    "Note I have modified the problem to now have pulse width $\\tau$ to better match the homework problem where $\\tau$ is a variable.\n",
    "\n",
    "The energy spectral density is\n",
    "\\begin{align}\n",
    "    G(f) = 4\\, \\text{sinc}^2(f\\tau)\\cdot\\sin^2(\\pi f\\tau)\n",
    "\\end{align}\n",
    "\n",
    "I convenient way to numerically integrate $G(f)$ is using simple reactangular partitions, but making sure that $\\Delta f$ is small relative to the changes in $G(f)$. Since you do not know what the value of $\\tau$ you consider a *normalized frequency* variable $f_n = f\\tau$ in the analysis. The rest of the steps are:\n",
    "\n",
    "1. Sweep $G(f_n)$ using an array `fn` running from zero to $f_n$ large enough to insure that $G(f_n)$ is very small relative to it largest value. In Python this is just filling an array, `Gn` with the functional values.\n",
    "2. Form a new array which contains the cumulative sum of the values in `Gn`, say `Gn_cumsum = cumsum(Gn)`. Aso form the sum of the array values, i.e., `Gn_tot = sum(Gn)`\n",
    "3. Plot the ratio of ``Gn_cumsum/Gn_sum` versus `fn`. The curve should start at zero and climb to one as $f_n$ becomes large. The value of $f_n$ where the curve crosses through 0.9 is the 90% containment bandwidth.\n",
    "\n",
    "**Note:** You might notice that $\\Delta f$, which is needed in the rectangular integration formula was never used. Why? In the calculation of the cumulative sum and the calculation of the total, both should include $\\Delta f$, hence in the ration the values cancel out. Nice!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fn = arange(0,10,.001)\n",
    "Gn = 4*sinc(fn)**2 * sin(pi*fn)**2\n",
    "Gn_cumsum = cumsum(Gn)\n",
    "Gn_tot = sum(Gn)\n",
    "plot(fn,Gn_cumsum/Gn_tot)\n",
    "grid()\n",
    "xlabel('Normalized Frequency $f\\tau$')\n",
    "ylabel('Fractional Power Containment');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fn_idx = np.nonzero(np.ravel(abs(Gn_cumsum/Gn_tot - 0.9)< 0.0005))[0]\n",
    "fn_idx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('The normalized 90 percent containment bandwidth is %2.2f Hz-s.' \\\n",
    "      % fn[1448])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Filter Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To facilitate the performance analysis of both discrete-time and continuous-time filters, the functions `freqz_resp()` and `freqs_resp()` are available (definitions below, respectively). With these functions you can quickly move from *z*-domain or *s*-domain rational system function coefficients to visualization of the filter frequency response\n",
    "* Magnitude\n",
    "* Magnitude in dB\n",
    "* Phase in radians\n",
    "* Group delay in samples or seconds (digital filter)\n",
    "* Group delay in seconds (analog filter)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def freqz_resp(b,a=[1],mode = 'dB',fs=1.0,Npts = 1024,fsize=(6,4)):\n",
    "    \"\"\"\n",
    "    A method for displaying digital filter frequency response magnitude,\n",
    "    phase, and group delay. A plot is produced using matplotlib\n",
    "\n",
    "    freq_resp(self,mode = 'dB',Npts = 1024)\n",
    "\n",
    "    A method for displaying the filter frequency response magnitude,\n",
    "    phase, and group delay. A plot is produced using matplotlib\n",
    "\n",
    "    freqs_resp(b,a=[1],Dmin=1,Dmax=5,mode = 'dB',Npts = 1024,fsize=(6,4))\n",
    "\n",
    "        b = ndarray of numerator coefficients\n",
    "        a = ndarray of denominator coefficents\n",
    "     Dmin = start frequency as 10**Dmin\n",
    "     Dmax = stop frequency as 10**Dmax\n",
    "     mode = display mode: 'dB' magnitude, 'phase' in radians, or \n",
    "            'groupdelay_s' in samples and 'groupdelay_t' in sec, \n",
    "            all versus frequency in Hz\n",
    "     Npts = number of points to plot; defult is 1024\n",
    "    fsize = figure size; defult is (6,4) inches\n",
    "    \n",
    "    Mark Wickert, January 2015\n",
    "    \"\"\"\n",
    "    f = np.arange(0,Npts)/(2.0*Npts)\n",
    "    w,H = signal.freqz(b,a,2*np.pi*f)\n",
    "    plt.figure(figsize=fsize)\n",
    "    if mode.lower() == 'db':\n",
    "        plt.plot(f*fs,20*np.log10(np.abs(H)))\n",
    "        plt.xlabel('Frequency (Hz)')\n",
    "        plt.ylabel('Gain (dB)')\n",
    "        plt.title('Frequency Response - Magnitude')\n",
    "\n",
    "    elif mode.lower() == 'phase':\n",
    "        plt.plot(f*fs,np.angle(H))\n",
    "        plt.xlabel('Frequency (Hz)')\n",
    "        plt.ylabel('Phase (rad)')\n",
    "        plt.title('Frequency Response - Phase')\n",
    "\n",
    "    elif (mode.lower() == 'groupdelay_s') or (mode.lower() == 'groupdelay_t'):\n",
    "        \"\"\"\n",
    "        Notes\n",
    "        -----\n",
    "\n",
    "        Since this calculation involves finding the derivative of the\n",
    "        phase response, care must be taken at phase wrapping points \n",
    "        and when the phase jumps by +/-pi, which occurs when the \n",
    "        amplitude response changes sign. Since the amplitude response\n",
    "        is zero when the sign changes, the jumps do not alter the group \n",
    "        delay results.\n",
    "        \"\"\"\n",
    "        theta = np.unwrap(np.angle(H))\n",
    "        # Since theta for an FIR filter is likely to have many pi phase\n",
    "        # jumps too, we unwrap a second time 2*theta and divide by 2\n",
    "        theta2 = np.unwrap(2*theta)/2.\n",
    "        theta_dif = np.diff(theta2)\n",
    "        f_diff = np.diff(f)\n",
    "        Tg = -np.diff(theta2)/np.diff(w)\n",
    "        max_Tg = np.max(Tg)\n",
    "        #print(max_Tg)\n",
    "        if mode.lower() == 'groupdelay_t':\n",
    "            max_Tg /= fs\n",
    "            plt.plot(f[:-1]*fs,Tg/fs)\n",
    "            plt.ylim([0,1.2*max_Tg])\n",
    "        else:\n",
    "            plt.plot(f[:-1]*fs,Tg)\n",
    "            plt.ylim([0,1.2*max_Tg])\n",
    "        plt.xlabel('Frequency (Hz)')\n",
    "        if mode.lower() == 'groupdelay_t':\n",
    "            plt.ylabel('Group Delay (s)')\n",
    "        else:\n",
    "            plt.ylabel('Group Delay (samples)')\n",
    "        plt.title('Frequency Response - Group Delay')\n",
    "    else:\n",
    "        s1 = 'Error, mode must be \"dB\", \"phase, '\n",
    "        s2 = '\"groupdelay_s\", or \"groupdelay_t\"'\n",
    "        print(s1 + s2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def freqs_resp(b,a=[1],Dmin=1,Dmax=5,mode = 'dB',Npts = 1024,fsize=(6,4)):\n",
    "    \"\"\"\n",
    "    A method for displaying analog filter frequency response magnitude,\n",
    "    phase, and group delay. A plot is produced using matplotlib\n",
    "\n",
    "    freqs_resp(b,a=[1],Dmin=1,Dmax=5,mode='dB',Npts=1024,fsize=(6,4))\n",
    "\n",
    "        b = ndarray of numerator coefficients\n",
    "        a = ndarray of denominator coefficents\n",
    "     Dmin = start frequency as 10**Dmin\n",
    "     Dmax = stop frequency as 10**Dmax\n",
    "     mode = display mode: 'dB' magnitude, 'phase' in radians, or \n",
    "            'groupdelay', all versus log frequency in Hz\n",
    "     Npts = number of points to plot; defult is 1024\n",
    "    fsize = figure size; defult is (6,4) inches\n",
    "    \n",
    "    Mark Wickert, January 2015\n",
    "    \"\"\"\n",
    "    f = np.logspace(Dmin,Dmax,Npts)\n",
    "    w,H = signal.freqs(b,a,2*np.pi*f)\n",
    "    plt.figure(figsize=fsize)\n",
    "    if mode.lower() == 'db':\n",
    "        plt.semilogx(f,20*np.log10(np.abs(H)))\n",
    "        plt.xlabel('Frequency (Hz)')\n",
    "        plt.ylabel('Gain (dB)')\n",
    "        plt.title('Frequency Response - Magnitude')\n",
    "\n",
    "    elif mode.lower() == 'phase':\n",
    "        plt.semilogx(f,np.angle(H))\n",
    "        plt.xlabel('Frequency (Hz)')\n",
    "        plt.ylabel('Phase (rad)')\n",
    "        plt.title('Frequency Response - Phase')\n",
    "\n",
    "    elif mode.lower() == 'groupdelay':\n",
    "        \"\"\"\n",
    "        Notes\n",
    "        -----\n",
    "\n",
    "        See freqz_resp() for calculation details.\n",
    "        \"\"\"\n",
    "        theta = np.unwrap(np.angle(H))\n",
    "        # Since theta for an FIR filter is likely to have many pi phase\n",
    "        # jumps too, we unwrap a second time 2*theta and divide by 2\n",
    "        theta2 = np.unwrap(2*theta)/2.\n",
    "        theta_dif = np.diff(theta2)\n",
    "        f_diff = np.diff(f)\n",
    "        Tg = -np.diff(theta2)/np.diff(w)\n",
    "        max_Tg = np.max(Tg)\n",
    "        #print(max_Tg)\n",
    "        plt.semilogx(f[:-1],Tg)\n",
    "        plt.ylim([0,1.2*max_Tg])\n",
    "        plt.xlabel('Frequency (Hz)')\n",
    "        plt.ylabel('Group Delay (s)')\n",
    "        plt.title('Frequency Response - Group Delay')\n",
    "    else:\n",
    "        print('Error, mode must be \"dB\" or \"phase or \"groupdelay\"')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: Discrete-Time Chebyshev Type I Bandpass Filter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sk_dsp_comm.iir_design_helper as iird\n",
    "import sk_dsp_comm.fir_design_helper as fird"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b1,a1,sos1 = iird.IIR_bpf(200,250,300,350,0.1,60.0,1000,'butter')\n",
    "b2,a2,sos2 = iird.IIR_bpf(200,250,300,350,0.1,60.0,1000,'cheby1')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "figure()\n",
    "iird.freqz_resp_cas_list([sos1,sos2],'dB',1000)\n",
    "ylim([-70,0])\n",
    "grid();\n",
    "figure()\n",
    "iird.freqz_resp_cas_list([sos1,sos2],'groupdelay_t',1000)\n",
    "grid();\n",
    "figure()\n",
    "iird.sos_zplane(sos2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b,a = signal.cheby1(5,.1,2*array([250,300])/1000,btype='bandpass')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "freqz_resp(b,a,mode='dB',fs=1000,fsize=(6,2))\n",
    "grid()\n",
    "ylim([-80,5]);\n",
    "xlim([100,400]);\n",
    "freqz_resp(b,a,mode='groupdelay_s',fs=1000,fsize=(6,2))\n",
    "grid()\n",
    "xlim([100,400]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: Continuous-Time Bessel Bandpass Filter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bc,ac = signal.bessel(7,2*pi*array([10.0,50.0])*1e6,btype='bandpass',analog=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "freqs_resp(bc,ac,6,9,mode='dB',fsize=(6,2))\n",
    "grid()\n",
    "ylim([-80,5]);\n",
    "freqs_resp(bc,ac,6,9,mode='groupdelay',fsize=(6,2))\n",
    "grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Second-Order Butterworth Lowpass Response"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider a 3rd-order analog Butterworth is the $s$-domain having transfer function $H(s)$. Using the `scipy.signal` function `butter()` we find the coefficients to the rational transfer function of the form:\n",
    "\\begin{align}\n",
    "    H(s) = \\frac{\\sum_{n=0}^M b_n s^n}{\\sum_{n=0}^N a_n s^n}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b3,a3 = signal.butter(3,2*pi*1,analog=True)\n",
    "freqs_resp(b3,a3,-1,2,mode='dB',fsize=(6,2))\n",
    "grid()\n",
    "ylim([-80,5]);\n",
    "freqs_resp(b3,a3,-1,2,mode='groupdelay',fsize=(6,2))\n",
    "grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Obtaining the Step Response via Simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Time domain simulation of continuous time system can be performed using the `signal.lsim()` function. You have to make sure the time step is sufficiently small relative to the filter bandwidth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = arange(0,2,.0001)\n",
    "xs = ss.step(t)\n",
    "tout,ys,x_state = signal.lsim((b3,a3),xs,t)\n",
    "plot(t,ys)\n",
    "title(r'Third-Order Butterworth Step Response for $f_3 = 1$ Hz')\n",
    "ylabel(r'Ste Response')\n",
    "xlabel(r'Time (s)')\n",
    "grid();"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}